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Pattern formation, i.e., the generation of an inhomogeneous spatial activity distribution 
in a dynamical system with translation invariant structure, is a well-studied phenomenon 
in neuronal network dynamics, specifically in neural field models. These are population 
models to describe the spatio-temporal dynamics of large groups of neurons in terms of 
macroscopic variables such as population firing rates. Though neural field models are often 
deduced from and equipped with biophysically meaningful properties, a direct mapping to 
simulations of individual spiking neuron populations is rarely considered. Neurons have 
a distinct identity defined by their action on their postsynaptic targets. In its simplest 
form they act either excitatorily or inhibitorily. When the distribution of neuron identities is 
assumed to be periodic, pattern formation can be observed, given the coupling strength 
is supracritical, i.e., larger than a critical weight. We find that this critical weight is strongly 
dependent on the characteristics of the neuronal input, i.e., depends on whether neurons 
are mean- or fluctuation driven, and different limits in linearizing the full non-linear system 
apply in order to assess stability. In particular, if neurons are mean-driven, the linearization 
has a very simple form and becomes independent of both the fixed point firing rate and 
the variance of the input current, while in the very strongly fluctuation-driven regime 
the fixed point rate, as well as the input mean and variance are important parameters 
in the determination of the critical weight. We demonstrate that interestingly even in 
"intermediate" regimes, when the system is technically fluctuation-driven, the simple 
linearization neglecting the variance of the input can yield the better prediction of the 
critical coupling strength. We moreover analyze the effects of structural randomness by 
rewiring individual synapses or redistributing weights, as well as coarse-graining on the 
formation of inhomogeneous activity patterns. 

Keywords: pattern formation, spiking neurons, linear model, mean-driven, fluctuation driven, ring networks, 
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1. INTRODUCTION 

Understanding the dynamics of neuronal networks and its depen- 
dence on connection topology, weight distribution or input 
statistics is essential to understand network function. One very 
successful field in uncovering structure-dynamics relationships 
are so-called neural field models (see e.g., Beurle, 1956; Wilson 
and Cowan, 1972; Amari, 1977; Ermentrout and Cowan, 1979a; 
Ben-Yishai et al., 1995; Ermentrout, 1998; Coombes, 2005). In 
these models the positions of neurons in space are substituted by 
densities and their respective coupling is described by coupling 
kernels depending on pairwise spatial distance. Given certain 
symmetries such as translation invariance, and homogeneous 
input the resulting non-linear dynamics is often reduced to a 
low-dimensional system that can be studied analytically. Such sys- 
tems can produce various spatio-temporal phenomena, such as 



traveling waves, activity bumps and formation of periodic pat- 
terns, which can be linked to activity in biological systems, e.g., 
visual hallucinations, feature selectivity, short term memory, or 
EEG rhythms (see e.g., Ermentrout, 1998; Coombes, 2005 and 
references therein). 

However, when thinking of actual neuronal tissue, the symme- 
try requirements and continuity assumptions apply on a rather 
macroscopic scale, when neuron densities are high and hetero- 
geneities are negligible. Moreover, such rate-based models cannot 
fully resolve the statistics of synaptic input currents, but it is 
well-known that neurons and neuronal populations react quite 
strongly to changes in the statistics of incoming currents in 
terms of mean, auto- and cross-covariance (see e.g., Mainen and 
Sejnowsky, 1996; Silberberg et al., 2004; Fourcaud-Trocme and 
Brunei, 2005; De la Rocha et al, 2007; Boucsein et al, 2009; 
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Tchumatchenko and Wolf, 2011). Finally, it is often unclear how 
to quantitatively interpret the parameters of abstract field models, 
especially when compared to spiking network simulations. 

Along these lines, Usher et al. (1995) demonstrated 
numerically how periodic pattern-formation takes place in 
a two-dimensional toroidal network of excitatory and inhibitory 
pulse-coupled integrate-and-fire oscillators. The bifurcation 
occurs when the relative surround inhibition reaches a certain 
threshold that is determined from the Fourier-transform of the 
effective coupling matrix (dispersion relation). In their study, 
Usher et al. (1995) moreover observed that the dynamics of the 
emerging pattern depends on the average external input current, 
with slowly diffusing hexagonally arranged activity bumps for 
low input current, over stationary bumps for intermediate input 
strength, to coherently moving bumps for strong external input. 

In a later study, Bressloff and Coombes (1998, 2000) showed 
rigorously how periodic orbits in networks of pulse-coupled 
inhibitory and excitatory neurons with Mexican-hat topology 
indeed undergo a discrete Turing-Hopf-bifurcation, if the system 
is coupled strongly enough. 

Roxin et al. (2005) studied a neuronal field model with partic- 
ular emphasis on the role of delays. Similar to earlier studies of 
neuronal field models, (see e.g., Ben-Yishai et al., 1995; Compte 
et al, 2000; Shriki et al., 2003) neurons were distributed along 
a ring and coupling was distance-dependent. In particular, the 
coupling kernel was of the form J(\x — y\) = J 0 + /icos(x — y), 
where x and y are neuron positions in space and Jo and J\ are cou- 
pling strength coefficients. Depending on the choice of Jo, Ji and 
the delay the system can show a plethora of activity states such 
as traveling waves, stationary or oscillatory bumps, and also ape- 
riodic behavior. However, this choice of interaction between two 
neurons is rather qualitative. It is usually symmetric and changes 
sign with distance, properties that are not in line with the biology 
of individual neurons and synapses. 

Thus, to compare the field-model results to simulations of net- 
works of excitatory and inhibitory spiking neurons, Roxin et al. 
(2005) set up a model where the same number of excitatory and 
inhibitory neurons were uniformly distributed across two rings 
and were coupled by sinusoidally modulated connection prob- 
ability. This set-up resulted in an effectively one-dimensional 
model where at each position excitation and inhibition combine 
to a net-coupling comparable to the field model, and spatio- 
temporal patterns qualitatively matched the field-model predic- 
tions. The parameters in the spiking model were chosen with 
regard to biophysically plausible values, and the non-linearity of 
the current-to-rate transduction in the field model was chosen as 
threshold linear in an ad-hoc manner, however a direct quantita- 
tive mapping between neuronal network and field model was not 
tried. 

Here, we investigate the dynamics of identical leaky integrate- 
and-fire neurons that are arranged on ring and grid networks 
where both inhibitory and excitatory neurons are uniquely 
assigned positions in space rather than representing a density. 
No two neurons can be at the same position, and this is thus 
a step toward spatially embedded networks of spiking neu- 
rons instead of neuronal populations with possibly ambiguous 
coupling weights. Instead, individual neurons have a distinctive 



identity in that they act either excitatorily or inhibitorily on all 
their post-synaptic targets (Dale's principle) and the coupling 
between any two neurons is usually not symmetrical. This is 
important since neglecting the identity of a neuron can lead to 
dramatically different dynamics (Kriener et al., 2008, 2009). If the 
system is translation-invariant, a Turing-instability occurs for a 
critical synaptic coupling strength J c , such that there is a bista- 
bility between the spatially homogeneous firing rate distribution 
and an inhomogeneous periodic spatial pattern, similar to what 
was described before (Usher et al., 1995; Bressloff and Coombes, 
1998, 2000). 

We find that the value of J c depends strongly on the statistics 
of the neuronal input: if neurons are mean-driven, i.e., receive 
suprathreshold input, the instability occurs for much smaller cou- 
pling strength than when neurons are in the strongly fluctuation- 
driven regime, where the mean input is subthreshold and spikes 
are evoked by spontaneous suprathreshold fluctuations. The latter 
is the key mechanism in explaining asynchronous-irregular firing 
in balanced neuronal networks. 

In assessing stability of an invariant network state to perturba- 
tions the usual modus operandi is to linearize around the fixed 
point (see e.g., Ermentrout and Cowan, 1979a,b, 1980; Usher 
et al, 1995; Bressloff and Coombes, 1998, 2000; Coombes, 2005), 
implying crucially the derivative of the neuronal input-output 
function. In the fluctuation-driven regime the input-output func- 
tion depends on both the mean and variance of the input current 
and thus both respective derivatives are expected to influence the 
quantitative prediction of linear fixed point stability (Amit and 
Brunei, 1997; Tetzlaff et al, 2012; Helias et al, 2013). Indeed, in 
the strongly fluctuation-driven regime, such that the standard- 
deviation of the input current is much larger than the average 
distance to firing threshold, the correct critical coupling strength 
is derived by linearization with respect to both mean and vari- 
ance around the fixed point. In the mean-driven regime, on 
the other hand, the linearization becomes independent on the 
input current variance, and the derivative with respect to the 
mean is constant over wide ranges of working points implying 
independence on the exact location of the rate fixed point. In 
this regime pattern formation typically occurs for considerably 
smaller coupling strength than in the fluctuation-driven regime. 
Interestingly, if the input current is in an intermediate regime 
with both moderate subthreshold input mean, yet considerable 
variance the linearization independent of the input variance can 
give a better estimate than the one taking it into account. We will 
see that usually, the true critical coupling strength will lie inbe- 
tween both estimates, with J c generally closer to the mean-driven 
regime prediction. 

Moreover, in order to relate our ring model to the earlier find- 
ings of pattern formation in the classical Ermentrout-Cowan net- 
works (Ermentrout and Cowan, 1979a,b, 1980) we also introduce 
a coarse-grained network that maps to a two-dimensional system 
which is formally identical to Ermentrout-Cowan networks, but 
lacks some of the details of the full, effectively five-dimensional 
ring network considered here. 

Finally, we study the role of robustness of pattern formation to 
randomness: first, structural randomness is caused by the intro- 
duction of short-cuts which move the ring topology into the 
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so-called small-world regime (Watts and Strogatz, 1998); sec- 
ondly, we also discuss the effect of randomness in the weight 
distribution introduced by lifting the constraint that neurons are 
either fully excitatory in their action or fully inhibitory (Dale's 
principle). 

2. MATERIALS AND METHODS 
2.1. NEURON AND NETWORK MODEL 

We study ring networks of N leaky integrate-and-fire (LIF) neu- 
rons with current-based synapses. Ne = fW, p € [0, 1], of all 
neurons are excitatory, the residual Ni = N — Ne neurons are 
inhibitory. The membrane potential Vj(f) of neuron i is governed 
by the differential equation 

dVi(t) 

x m — p = -Vi(t) + RI S j(t -d) + W X)i (f) , (1) 
at 

where the membrane potential Vj(f) is reset to V res whenever 
it reaches the threshold potential Vthr and a spike is emitted. 
The neuron is then refractory for a time T re f. t m denotes the 
membrane time constant, RI s ,i(t — d) = x m Ylf= [ WgSj(t — 
d), where / s is the input current received from the local net- 
work, and s;(t) = J^k — f i,*) ^ s tne spike train produced by 
some neuron j. d is the transmission delay, and W is the coupling 
matrix, with entries that are either zero (no synapse present), /e = 
/ (excitatory synapse) or ]\ = —gj (inhibitory synapse). To keep 
spiking activity from entering a high rate state, we assume that 
the recurrent input I s is net inhibitory by choosing g > Ne/Ni. 

I x denotes the external input current that is modeled as 
stationary Poisson noise injected via current-based synapses 
of strength / x (for details, see below), and R is the membrane 
resistance. Such leaky integrate-and-fire neurons were extensively 
studied in various input (see e.g., Amit and Tsodyks, 1991; Amit 
and Brunei, 1997; Lindner, 2004; Fourcaud-Trocme and Brunei, 
2005; De la Rocha et al, 2007; Vilela and Lindner, 2009; Helias 
et al., 2010) and network settings, especially in the context of 
balanced random networks (Brunei, 2000; Tetzlaff et al, 2012). 
Here, we consider ring networks of LIF neurons, such that each 
fifth neuron is inhibitory and the coupling matrix W is given by 
(cf. Figure 3 A) 



two-dimensional model (see e.g., Ben-Yishai et al, 1995; Compte 
et al., 2000; Shriki et al., 2003; Roxin et al, 2005). As we will 
demonstrate in section 3.2, the network studied here leads to an 
effectively five-dimensional model. In section 3.4 we will discuss 
a reduction to a two-dimensional model that however already 
deviates noticeably in its details from the full five-dimensional 
network. 

2.2. CONSTANT EXTERNAL INPUT IN THE MEAN-DRIVEN AND 
INTERMEDIATE REGIME 

In the case of the mean-driven (|i[i^fj > Vthr — ^res) and inter- 
mediate regimes, which will be discussed later, the external 
current is purely excitatory, such that RI x (t) = t m / x s x (f) with 
/ x > 0 and rate vx = E[s x (f)L where E[.] denotes the expec- 
tation value. To quantify the effective strength of the external 
input, we express vx = n9// x T m , i.e., the parameter r\ gives the 
strength of the external drive in terms of the excitatory rate nec- 
essary to reach threshold. Assuming stationary spiking activity 
with rate v 0 for both the excitatory and inhibitory population, 
i.e., ve = E[s£(f)] = v 0 and vi = E[si(f)] = v 0 respectively, the 
mean and variance of the total input RI = R(I S + J x ) are thus 
given by 

H [R (I s + 7 X )] = £ X m VjJj = KT m v 0 /(P - g{l - P)) 
;€{E,I,X| 

+ tmVx/x (3) 

a 2 [R (I s +k)]= Tm V, 2 = ^mV 0 / 2 (P + g 2 (l - P)) 

je{E,I,X] 

+ tmVx/ x - 

For stronger local coupling strength / the negative mean recur- 
rent contribution u.[W s ] to the total input is stronger, while the 
variance a 2 [RI s ] increases. For the networks we consider here, 
for small / the network activity is thus typically driven by the 
mean [i[RI] of the total input, while for larger / the variance 
a 2 \RT\ has a notable impact. Moreover, the network firing rate 
v 0 depends on /. 



W„ 



if j excitatory, \ i ■ 
if j inhibitory, \i - 
otherwise 



-;'|modN < k/2 and i ^ j 
- j\ modN < k/2 and i ^ j ■ 

(2) 



Thus, each neuron is connected to its k nearest neighbors, where 
we assume k = 0.1N. We note that the network layout chosen 
here is motivated by an actual embedding of individual neu- 
rons into some space, in contrast to earlier studies of pattern 
formation in spiking neuron networks where inhibitory and exci- 
tatory neurons are often distributed in equal numbers along two 
identical rings, and where the main difference lies in the spa- 
tial interaction ranges. These systems are usually motivated by 
population-density descriptions and can often be reduced to a 



2.3. EXTERNAL INPUT ADJUSTED TO ENSURE CONSTANT TOTAL 
CURRENT IN THE STRONGLY FLUCTUATION-DRIVEN REGIME 

The strongly fluctuation-driven regime is characterized by 
subthreshold mean total input, |x[_RI] < V t h r — V^s, but high 
variance a [RI\, such that spikes are initiated by transient 
input fluctuations. To guarantee that we are in the high 
variance regime for all /, we inject excitatory and inhibitory 
external currents with rates ve x , vi x , such that the total 
input current mean and variance, i.e., \i[RI] and a [RT\, 
stay the same when varying the local coupling strength /. 
This is achieved by adapting u.[-R-f x ] = \i[RI] — |x[W s ] and 
c 2 [^x] = o 2 \RT\ — 0 2 [_RI s ] as a function of ve x , vl x . The 
external current injected into each local neuron by current- 
based synapses of strengths /e x = / x and /ix = —gf x is thus 
given by 
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RI x (t) = W x (s Ex (f) -£ST x (f)) 

1 /CT 2 [£/ X ] 



(4a) 



with v Ex := E[s Ex (t)] 



and vix := E[s Ix (0] 



Wx(l +g) V /x 



(4b) 



1 



a 2 [RI x 



(4c) 



Because the variance of the recurrent input cr 2 [.RJ s ] increases 
quickly with increasing /, the total variance o 2 [Rr\- in 
order to keep it fixed for all /- usually needs to be cho- 
sen quite large in order to avoid non-sensical rates vt x < 0. 
This moreover implies that this approach of choosing 



the external input is not useful for the mean-driven 
scenario. 

For all neurons we set 9 = Vthr — V les = 20 mV, t m = 20 ms, 
and -R = 80 MQ,. The excitatory synaptic coupling strengths /, / x , 
the relative strength of inhibition g, as well as x re f and d will be 
specified individually. The model and parameters are listed in 
Table Al in Appendix Al. All simulations were carried out with 
NEST (Gewaltig and Diesmann, 2007). 

3. RESULTS 

For small coupling strength / and large mean current input 
|x[_Rf] = E[R(I X + I s )] the spiking dynamics of the ring net- 
works is characterized by locally clustered activity on the spa- 
tial scale of the footprint k, cf. Figures 1A,G (Kriener et al., 
2009). For increasing /, however, a clear spatial pattern emerges, 
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FIGURE 1 | Pattern formation in a ring network of 2500 neurons 
with g = 6 and varying absolute coupling strength parameter J. 
(A-C) Show the mean-driven case, where the external input is given by 
Poissonian input spikes with rate v x = 106/J x x m and amplitude 
J x = 0.1 mV, resulting in an excitatory input current of amplitude 
/ x = 2500pA. The critical coupling strength J c for the spatially 
homogeneous rate distribution to become bistable with a spatially 
modulated pattern is around 0.5 mV. (D-F) Show the same network in 
the fluctuation driven regime, where mean and variance of the external 
input current l x were chosen such that mean and standard deviation of 
the total input current l=l x + l s were fixed at yi[RI] = 5m\/ and 
a[RI] = 60 mV, cf. Equation (4). (G-l) Show the corresponding plots for 
an intermediate input regime, where the external input current was 



Poissonian with rate 3.59/J x x m and amplitude J, = 0.1 mV, resulting in 
an input current of amplitude / x = 875pA. Here, refractory time constant 
T re f and delay d were chosen to be 0.1 ms to minimize loss of input 
current and avoid pronounced delay oscillations in the mean-driven case 
(Helias et al., 2013). As we will demonstrate in section 3.2, eigenvalue 
analysis predicts an instability of an eigenmode with non-zero 
wavenumber A c > 0 (here, A c = 13) for coupling JJT d »;0.5mV in the 
mean-driven case, and j£ d «i0.9mV in the fluctuation-driven case. We 
see that already for J < J c there can be deviations from the spatially 
homogeneous mode A = 0 due to the fluctuating nature of the activity 
and selective amplification (Dayan and Abbott, 2001), especially in the 
fluctuation-driven case. For larger J the pattern becomes more 
prominent and finally some neurons fall completely silent (C,F,I). 
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cf. Figures 1C,F. Since all neurons receive by construction sta- 
tistically identical input, all neurons should fire at the same 
rate. But there is a distinct critical coupling strength J c beyond 
which the homogeneous distribution of firing rates becomes 
unstable and the system enters a spatially inhomogeneous state 
after sufficient perturbation (Bressloff and Coombes, 1998). We 
note that the spatio-temporal spiking activity can be quite rich 
in networks with translation-invariant symmetry, as discussed 
for example in (Ben-Yishai et al, 1995; Usher et al., 1995; 
Bressloff and Coombes, 1998, 2000; Shriki et al, 2003; Roxin 
et al, 2005; Marti and Rinzel, 2013). Here, we do not system- 
atically analyze the temporal aspects of the activity in terms 
of traveling waves or oscillations, but focus on the onset of 
periodic pattern formation in space as a function of coupling 
strength. 

What determines this critical coupling strength J c where the 
bistability of the homogeneous spatial rate distribution and the 
spatially modulated distribution occurs? We find that J c depends 
strongly on the input regime the neurons operate in: if neurons 
are driven predominantly by the mean of the input current, i.e., 
if \l[RI] > 9, we find that J c is much smaller than if neurons 
are strongly fluctuation-driven, i.e., [i[RI] < 9 with large enough 
variance of the input current to occasionally drive the membrane 
potential to threshold, cf. Figure 1. 

In the following we present two different linear rate models 
of the integrate-and-fire dynamics that explain and predict the 
respective occurrence of the pattern formation. 



respective linearization in the limits of subthreshold mean input 
with fluctuation-driven spiking in section 3.1.1, and of domi- 
nant suprathreshold mean total input current |x[_R(/ x + 7 S )] > 9 
in section 3.1.2. In the fluctuation-driven regime a quantitative 
approximation is obtained from the modulation of both mean 
and fluctuations of the synaptic input to the neurons, while in the 
mean-driven regime we recover the intuitive result that the sys- 
tem is governed by the noise-free input-output relation given by 
the /-/-curve of the individual neuron. In general, the full spiking 
dynamics will lie somewhere in between these two limiting cases, 
the case we call "intermediate." 

3. 1. 1. Linearization of the input-output relation in the 
fluctuation-driven regime 

We first discuss the case when the system is fluctuation-driven, 
i.e., the average input current is subthreshold (|x[_RI] < 9) and 
spikes are induced by membrane-potential fluctuations. 

To assess linear stability in this regime we need to derive 
a linear model that self-consistently describes the dynamics of 
the fluctuations M,(f) = (s,(f) — v 0 ) = v,(r) — v 0 of the activity of 
neuron i around its mean firing rate or working point v 0 . Treating 
the response of integrate-and-fire neurons by first order pertur- 
bation theory is a well-established method first introduced in 
Abbott and van Vreeswijk (1993) for the perfect integrator and in 
Amit and Brunei (1997); Brunei and Hakim (1999) for the leaky 
integrate-and-fire neuron. 

The dynamics of u(t) can then be written in the form 



3.1. TWO LINEARIZATIONS FOR THE SELF-CONSISTENT RATE IN 
NETWORKS OF LIF-NEUR0NS 

The analysis of the stability of spatially homogeneous activity 
dynamics to spatial perturbations follows a simple concept: first, 
the stationary, homogeneous state is determined in a mean-field 
approximation. In particular, in the diffusion limit, i.e., under 
the assumption that the coupling strength / is small compared 
to the distance between reset and threshold 9 and that all neurons 
receive statistically identical input, the stationary firing rate v 0 of 
the neurons can be derived self-consistently as a function of mean 
and variance of the input current (a solution of the mean first- 
passage-time problem first derived by Siegert (Siegert, 1951), but 
also see e.g., Amit and Tsodyks, 1991; Brunei, 2000), such that 

v^ 1 = T re f + T m .y/tt / exp \x 2 ] (1 + erfix]) dx, (5) 

/ Vres - V-o L J 



u(f) ~ h(t) * [u(t) + x(t)] , (6) 

where h(t) is a linear filter and x(t) is a white noise imitating 
the spiking nature of the signal. In the Fourier-domain Equation 
(6) becomes the simple product L7(orj) = H(u>)[U(u>) + X(w)], 
where the capitalizations indicate the Fourier-transforms of h(t), 
u(t) and x(t). Pattern formation is typically a slow process and 
the transfer function H(w) of leaky integrate-and-fire neurons 
is dominated by low-pass behavior \ In order to assess stability, 
we therefore estimate H(co) in the low-frequency limit, i.e., 

H(0) = f Q °° h{t) dt. With = x m W I} and ^ = t m W| as 

well as = 2o 0 j^£r, we obtain (see also Amit and Brunei, 
1997; Tetzlaff et al, 2012; Helias et al, 2013 for details) the 
effective coupling matrix W of the system as the local derivative 
of the input-output-curve given by Equation (5): 



where |x 0 = £ie{E,i,x}70'«>' c m and o 2 0 = £ je |E I X| jfv io x m are 
the self-consistent input mean and variance. Because all neu- 
rons receive statistically identical input, in absence of symmetry- 
breaking the firing rate distribution is spatially homogeneous, 
i.e., v 0 j = v 0 for all i. To assess simple stability of the homoge- 
neous state to spatial perturbations, we need to perform linear 
perturbation analysis. The critical eigenvector of the resulting 
linear stability operator then determines the spatial pattern that 
develops. 

We will show that the exact nature of the appropriate lineariza- 
tion depends critically on the input current regime, and derive the 



WyiWij) = / hyit) dt=-^ (7a) 
Jo dSj 

= ^ Xm W lj+ -^-^Wl (7b) 



Note, that for exponential response kernels h(t) = x exp[— t/x] the response 
is structurally equivalent to a linear first order differential equation, though it 
has been demonstrated that in ring networks this response kernel can be sub- 
stantially different from exponential and even contain oscillatory components 

(Lindner and Schimansky-Geier, 2001; Pernice et al., 2012). 
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(v 0 t m ) 2 V^ (/(y e ,i) (l + -^ye,; 
o 0 ,i \ \ 2a 0 j 

( Wit 
-/Ov.i) 1 + ; — Yr.i 
\ 2a o j 



(7c) 



with 



f(y) = /(erfty] + 1), y e = — , and y r = ^—^.(8) 



The effective coupling matrix W(Wy) is thus structurally identi- 
cal to W given by Equation (2), with Wy = W(J) if j is excitatory, 
and Wy = W(—gJ) if j is inhibitory, yielding an effective rela- 
tive inhibitory strength g := \W(—gJ)/W(J)\. The zero frequency 
mode 1/(0) therefore fullfills the equation 



1/(0) = [I - W] 1 X(0) . 



(9) 



We expect linearity to break down if W has any eigenvalue with 
real part larger than unity. The critical coupling strength J^ d 
for which this is the case can be determined implicitly from 
Equation (7). 

3. 1.2. Linearization in the low input current noise limit 

Here, we discuss the linearization in the case that neurons are 
mean-driven, i.e., when the mean of the input current [i[RI] is 
suprathreshold, while the variance of the input is negligible [see 
also Amit and Tsodyks (1991) for analogous arguments]. In this 
case, the input-output-relation of the neuron becomes effectively 
linear in [i[RI] until it reaches saturation at v max = l/x re f (Amit 
and Tsodyks, 1991; Salinas and Sejnowsky, 2000). This can be seen 



by considering the input-output relation of an individual neuron 
in response to a constant current, i.e., the /-/-curve 



v(I) 



tref - TmLog 



RI 



(10) 



For RI 3> 6, i.e., Jj <SC 1 and negligible t re f, it holds (see 
Appendix A2 for details) 



v(I) 



RI 

T m 8 



1 

2x„ 



(11) 



This linear- threshold-type relationship v(I) = — ^— J fits 

that observed in simulations very well, cf. Figures 2A,B. 

To derive a rate equation for the full spiking network dynam- 
ics, we first define the instantaneous rate of neuron i as 

1 Mt/2 

Vj(t) = E A( [s,-(0] := lim At ^ 0 — / E[s t (t)]dt, (12) 

where E[.] is the average over realizations. 

Substituting the spike trains by their respective firing rates, the 
input to neuron i takes the form RI, = i m W,jVy + RI X . If we 
only consider networks with rates v,- > 0 and approximate the fir- 
ing rate of neuron i by the linear expression (11), we arrive at the 
linear equation 



wn 1 or/ x - e/2) . 



(13) 




[mV] J[mV] J[mV] 



FIGURE 2 | (A) Demonstrates that the input-output-reiation of the LIF neuron 
Equation 5 indeed gets linear in the strongly mean-driven regime. The light 
gray line shows the f-/-curve in the noise-less case [Equation (10)], dark gray 
corresponds to its linear approximation Equation (11), while the red curve is 
the corresponding self-consistent rate given by Equation (5) in the low-noise 
limit a[RI] -> 0. (B) Shows the output rates as a function of network coupling 
strength J as they are actually obtained from network simulations averaged 
over 10 trials (blue dotted curve) vs. the prediction from Equation (5) (gray 
left-pointing triangles), as well as the linear model prediction Equation (13) 
with E[l/(f)] = 0/2 (dark gray triangles). Also shown is the prediction from 
solving Equation (14) self-consistently with V = f_ VP(V) dV [with 
stationary membrane potential probability density function P(V) as derived in 
Brunei (2000); gray right-pointing triangles], as well as from using the 



estimated average membrane potentials E[V/(f)] obtained from simulations 
(black asterisks). The vertical red line represents the expected critical 
coupling strength J md for the mean-driven regime. The external drive is 
constant Poissonian noise (ri = 3.5), while the local network coupling 
strength J, and hence m-I^'s] and o[RI B ], vary. In this setting the system 
undergoes a regime change from the mean-driven to the fluctuation-driven 
scenario. This is demonstrated in (C), where the total mean and standard 
deviation of RI = R(/ x + 4) are shown as function of J in black and gray, 
respectively. The red dashed line corresponds to that in (B), while the black 
dashed line indicates the critical expected from Equation (9). The 
corresponding spike activity is shown for three exemplary cases in 
Figures 1G-I. Other parameters in (B,C) are N = 2500, k = 250, 8 = 20 mV, 
x m = 20 ms, x ref = 0.1 ms, and g = 6. 
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Indeed, if t re f is negligible, Equation (1) can be rewritten as 
(Kriener et al, 2008) 

dViit) 

x m —^ = -Vi(t)-x m Qsi(t) + x m 22 W tjS j(t - d) + RI Xii . 

(14) 

Equation (13) can then be derived analogously by tempo- 
ral averaging, assuming stationarity dE& t [V(t)]/dt = 0, i.e., 
EAtfVXf)] = const, and thus EAtfv(f)] = const = v 0 . For high 
rates in response to large |x[iU] we expect E[V(f)] ~ (Vthr — 
Vres)/2 = 6/2, as the voltage passes through the interval 
between reset and threshold with approximately constant velocity, 
leading to 

o = -e/2 - t m ev,(f) + im J2 w, j v j + (15) 

i 

which is identical to (13). We note, that the self-consistent solu- 
tion of Equation (15) yields the correct quantitative rates over 
a wide range of relative input magnitude, if instead of 6/2 
the actual E[V r ,(t)]-values measured in simulations are inserted. 
If all neurons are identical and receive statistically identical 
input, this mean value can be obtained from V,- E[ V,(f)] = V := 
/® VP(V) dV, where P(V) denotes the stationary membrane 
potential probability density function as e.g., derived in Brunei 
(2000). All output-rate predictions are compared to the out- 
come from simulations in Figure 2B. Aside from the prediction 
that assumes E[V(f)] = 8/2 (dark gray triangles) and underes- 
timates the true rate, all predictions fit the simulation results 
(blue dots) very well. In particular, the Siegert equation (5) 
coincides perfectly with the linearized rate Equation (14) if we 
assume E[V(f)] = V (light gray triangles), while the best fit is 
obtained with Equation (14) and the actual measured E[V,(f)] 
(black asterisks). 

The stability of the fixed point rate to small perturbations is 
determined by the largest real part of all eigenvalues of W/Q. If 
one eigenvalue X c has real part larger than unity, we expect the 
corresponding eigenvector v c to grow exponentially, only limited 
by the non-linearities due to rate-rectification for small rates, and 
the saturation because of neuronal refractoriness for high rates. 
This determines the critical coupling strength 



that in this limit is independent of the exact working point v 0 > 0. 
Even though for coupling strengths / > /™ d the rate instability 
can lead to pattern formation such that individual neurons fire 
at quite different rates, the population rate is still captured well 
by the firing rate predictions assuming homogeneous rates across 
neurons, see Figure 2B. 

Finally, we note that for the constant external drive scenario 
considered here the mean and variance of the total input cur- 
rent vary with varying /. In particular, for increasing / the system 
undergoes a cross-over from the mean-driven to the fluctuation- 
driven regime, see Figure 2C. 



3.2. EIGENSYSTEM OF DALE-CONFORM TRANSLATION-INVARIANT 
RING NETWORKS 

In this section we will analytically derive the eigenvalue spec- 
trum of the ring networks under consideration that yields the 
critical eigenvalue and thus the critical coupling strength J c with 
respect to the two linearizations Equation (9) and Equation (13). 
As described in section 2 we consider ring networks of size N 
with regular connectivity, in that each neuron is connected to 
its k nearest neighbors, k/2 on each side of the neuron but 
not to itself. Inhibitory neurons shall be distributed periodically 
across the network as illustrated in Figure 3A, where the ratio 
between excitation and inhibition is 4 : 1 (|3 = 0.8). This is in line 
with the observed frequencies of excitatory and inhibitory cells 




Wavenumber A=No/2tt£ Neuron i 



FIGURE 3 | (A) Sketch of the neuron layout on the ring topology. Each 
neuron is connected to its k nearest neighbors. A shift in neuron label by 
five yields the same ring as before and is formally expressed by the shift 
operator T t with 1 = 5. (B) Shows the coupling matrix W of the ring: black 
squares depict inhibitory, white squares excitatory, and gray squares zero 
coupling from / to / (W = 60, k = 30). (C) Eigenvalue spectrum of a ring 
matrix W/6 with N = 2500, k = 250, J = 1 mV, and g = 6. (D) The two 
eigenvectors belonging to the twice degenerated eigenvalue from (C) with 
largest real part l c . (E) The eigenvalues form five bands and are shown as a 
function of the wavenumber. The black curve depicts the band containing 
the critical (maximal) eigenvalue X c which here corresponds to wavenumber 
A c = 13, so the two crictical eigenvectors have 13 major peaks (cf. D,F). (F) 
The rate as a function of the neuron index in a simulation of N = 2500 
integrate-and-fire neurons with J = 1 mV and l x = 750 pA. Note that the 
spectra (shown in C) look slightly different in the fluctuation-driven case, 
because the absolute amplitude of the entries \ Wf\ j= \ and also the 
balance between positive and negative entries g j= g. The maximum of the 
dispersion relation and the eigenvectors are, however, the same. 
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J[mV] 

FIGURE 4 | Eigenvalue of W/d (red) and W (varying gray levels) with 
largest real part \ c in dependence of J and external drive n (increasing 
m = (3.5, 5, 10, 15, 20, 30, 40) from dark to light gray). The eigenvalues 
of l/V/6 do not depend on r| and thus the J c value such that X c = 1 
[indicated by red circle, cf. Equation (13)] holds uniquely for J™ d = 0.506. 
The gray lines show the dependence of \ c of W as obtained from Equation 
(7) with both contributions from dv/d\i., dv/do (solid lines) and with the 
dv/d\i contribution only (dashed, here for clarity we just plot the curve for 
r| = 3.5). The black triangle and the gray asterisk indicate the respective 
predictions of Jj. d = 1 .54 and j™ d - dv /* = 0.89 for t] = 3.5. For comparison, 
the blue line depicts the J-dependence in the input scenario with fixed total 
input current mean and variance, i.e., \i[RI ] = 5 mV and a[R/J = 60 mV, 
independent of changes in J. In that case, the self-consistent population 
firing rate is expected to be constant and thus Equation (7) is of the form 
W = c-\ Wjj + C2 W~, with constants Ci , C2- For large a and c-\ > C2 the 
curve is dominated by the linear weight dependence. For the parameters 
shown here, the expected Jj. d = 0.905 mV (cf. Figure 1). 



in cortex (Schtiz and Braitenberg, 1994). The ring can thus be 
divided into structurally identical elementary cells of four exci- 
tatory neurons and one inhibitory neuron. We choose k as an 
integer multiple of the size I of an elementary cell (here I = 
5, Kmod£ = 0), such that it moreover divides the total num- 
ber of neurons N (iVmodic = 0). The resulting coupling matrix 
is sketched in Figure 3B, where black squares depict inhibitory 
weights, white depict excitatory weights and gray denotes the 
lack of a connection. In this way all neurons receive the exact 
same amount of inhibition and excitation, and the translation 
invariant mode (1, 1, • • ■ , 1) T is an eigenvector of the coupling 
matrix with eigenvalue \i = k/(P — g(l — P)). To keep network 
activity balanced, inhibitory synapses are assumed to be g-times 
stronger, with g > 4, cf. Equation (2). Hence, the local network is 
inhibition dominated. 

To compute the eigensystem of the N-dimensional system we 
make use of its symmetry properties. The coupling matrix com- 
mutes with the unitary operator that shifts all neuron indices by 
multiples of I = 5. Hence, we can diagonalize both in the same 
eigenbasis and moreover reduce the problem to an effectively five- 
dimensional one as outlined in Appendix A3. Figure 3C shows 
the exemplary eigenvalue spectrum of in this case the rescaled 
coupling matrix W/Q of a network of size N = 2500 with absolute 
coupling strength / = 1 mV and relative strength of inhibition 
g = 6. Thus, the eigenvalue \ c of W/Q with largest real part will 
exceed unity, if the amplitude of the absolute coupling strength / 
exceeds 7 c md = l/Re[X c ], cf. Equation (16), here / f md « 0.5. 

The corresponding critical weight J { c d for the fluctuation- 
driven case is obtained by the identical eigenvalue decomposition 
of W as defined by Equation (7). Since W depends non-trivially 
on the self-consistent working point v 0 of the system, the critical 
coupling strength J™ is the solution of 

jf ■- J such that max,{Re[Xi(W(\><,, /))]] = 1 , (17) 

where X;, i e {1, . . . , N} are the N eigenvalues of W. If/ > Jf d/fd 
(depending on the regime), the spatially homogeneous rate distri- 
bution becomes unstable to perturbations in favor of the fastest 
growing eigenmode v c , which only depends on the symmetry 
properties of the matrix as long as g and g are larger than 
four, and it is hence the same irrespective of the input regime, 
cf. Figures 1C,F,I. The spatial frequency of the emerging spa- 
tial pattern (see Figure 3), i.e., the number of maxima of the 
rate distribution along the ring, is thus given by the wavenum- 
ber A (see Appendix A3) that corresponds to the respective 
X c := max,[Re[X,]]. We note, that the eigenvalue spectra in the 
fluctuation-driven case look slightly different from those in the 
mean-driven case (cf. Figure 3C) because of the in general dif- 
ferent scaling g / g between positive and negative entries in W. 
The maximum of the dispersion relation determining A c and the 
eigenvectors are, however, usually the same. 

Figure 4 shows the critical eigenvalue X c as a function of / 
for the linearizations Equation (9) (gray lines) and Equation 
(13) (red line) in the constant external drive regime, parameter- 
ized by T). The critical eigenvalue J c is given by the intersection 
with the horizontal line at one. As was pointed out in sec- 
tion 3.1.2, in the noiseless approximation Equation (13) the 



eigenvalues of the linear operator do not depend on the work- 
ing point, and thus the critical weight /™ d is unique for all x\ 
such that v 0 > 0 (here, critical /™ d = 0.506 indicated by red 
circle). 

For the linearization Equation (9), and taking into account 
both |x and a in Equation (7b) (solid lines), the criti- 
cal weight J { c d is strongly dependent on the value of r\ = 
{3.5, 5, 10, 15, 20, 30, 40} (increasing from dark to light gray), 
and never comes close to the prediction of the noiseless lineariza- 
tion /™ d . In particular, for larg er r| jf increases again. 

If we neglect the cr-dependence in Equation (7b), and just 
take into account dv 0 /d[i 0 = dvf A /dRI = l/t m 9 as given by 
Equation (11), Equation (7) is linear in Wy and we recover the 
linear noiseless approximation Equation (13) for Equation (9). 
We might thus expect that Equation (9) in general will give a J { c d - 
prediction close to /™ d if we neglect the cr-dependent quadratic 
term in Equation (7b). Although the prediction of this f c d \ dv °/ d v-° 
becomes smaller (cf. black triangle for J { c d = 1.54 mV vs. the gray 
asterisk for f^ dVi >/ d v-° = 0.89 mV for r| = 3.5) it never becomes 
as small as /™ d . This is because every change in [L„ will also 
influence a 0 via its impact on v 0 , and thus dv 0 ([i 0 , a 0 )/d[i 0 ^ 
dv™ d (RI)/dRI. 

We note that the observation /^ d > f^' dv °l d ^° [ s explained by 
the fact that W(W«)ldv 0 /rfu, 0 is linear in Wy, and thus the ratio 
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g = | W(-gJ)/W(J) I = g, while for the full Equation (7), W(W,j) 
is quadratic in Wij, yielding &g<g- Evaluating Equation (A20) 
in Appendix A3 shows that the critical eigenvalue \ c in both lin- 
earizations is a linearly increasing function of g, g, respectively, 
and J c is thus decreasing as a function of g, g. 

For comparison, Figure 4 also shows the /-dependence of the 
real part of the critical eigenvalue X c in the strongly fluctuation- 
driven regime ( |jl = 5mV, a[RI] = 60 mV, blue line), with 
the critical J c indicated by the crossing of the unity-line (blue 
triangle). 

3.3. LINEAR STABILITY IN DEPENDENCE OF THE INPUT CURRENT 
REGIME 

We will now compare the general performance of the lineariza- 
tions Equations (9) and (13) and in particular the resulting 
predictions of the critical coupling strength J c with the actual 
onset of pattern formation as observed in simulations. Figure 1 
demonstrates how the identical ring network structure undergoes 
pattern formation at very different values of network coupling 
strength / depending on the input current regime. In particular, 
pattern formation sets in later if the system is in a fluctuation- 
driven regime. This can be understood in the two limit-cases that 
were analyzed in sections 3.1.1 and 3.1.2. 

But what determines pattern formation onset in intermediate 
cases that comprise the most interesting and relevant cases? As we 
saw in section 3.1.2, the network will even undergo a crossover 
from the mean-driven to a more fluctuation-driven regime, if the 
external input is kept constant and not adapted to counterbal- 
ance the changes in recurrent input structure for changing /, see 
Figure 2C. Also, due to the spiking nature of the activity, even in 
a strongly mean-driven scenario there will always be a consider- 
able amount of variance which might change the onset of pattern 
formation with respect to /™ d . 

To give more insight into the actual onset of pattern formation, 
useful reduced measures are the variance and kurtosis of the dis- 
tribution of firing rates over neurons, cf. Figure 5. Figures 5A-C 
show the histogram of rates, and the variance and kurtosis of 
these histograms for various /-values in the mean-driven regime 
averaged over 10 trials each. Figures 5 D-F show the same for 
the fluctuation-driven case and Figures 5G-I for the intermediate 
case, all for the same parameters as in Figure 1. 

When the system is still well in the linear regime J < J c the 
distribution is approximately Gaussian with small variance and 
a kurtosis close to zero (cf. Figures 5A,D,G, light gray curves), 
while once the critical weight ] c is surpassed the variance increases 
strongly (cf. Figures 5B,E,H) and the kurtosis becomes negative 
(platycurtic), indicating the broadening of the distribution due 
to the inhomogeneous rate per neuron (cf. Figures 5C,F,G and 
5A,D,G, dark gray curves). The respective critical weights are 
indicated by vertical dashed lines in Figures 5B,E,H and C,F,I for 
visual guidance (red: / c md , black: / c fd , gray: f^ dv »/ d ^)_ 

The onset of pattern formation in the strongly fluctuation- 
driven regime at y~c is more shallow than in the mean-driven 
regime at /™ d . This might be due to deviations of the input spike 
statistics from the asynchronous-irregular activity assumption 
underlying the linear response derivation of /J d : the standard- 
deviation of the input is extremely large at a = 60 mV. This 



implies that the membrane-potential will be hyperpolarized for 
long times, while at other times it is highly depolarized, and 
several spikes are emitted in short succession. This typically 
leads to higher coefficients of variation, i.e., more irregular fir- 
ing than expected for Poisson spike statistics [see also Brunei 
(2000) and supplementary material section SI]. Also, even in the 
highly fluctuation-driven regime there may be residual correla- 
tions between spike trains. Lastly, the increased network coupling 
strength might lead to deviations from the Gaussian white noise 
approximation of input currents underlying Equation (5) and the 
linear response theory yielding Equation (7) (see supplementary 
material section SI for an analysis). How such deviations of the 
input statistics can indeed change the spike response behavior of 
LIF neurons was studied, e.g., in Moreno et al. (2002); Renart et al. 
(2007); Moreno-Bote et al. (2008); Helias et al. (2010). 

However, another contributing factor, and more likely the rea- 
son for the shallow transition at is the degeneracy of the 
maximal eigenvalue, see Figures 3D,E. In the fluctuation-driven 
regime the system is subject to strong perpetual perturbations 
that lead to a switching between the two dominant eigenvectors 
depicted in Figure 3D in the transition regime jf ±AJ. Indeed, 
as can be seen from Figure IE, the periodic pattern is already 
clearly distinguishable at J c d over several hundred milliseconds, 
but the average activity depicted in the histogram is still quite 
flat, yielding smaller variance and kurtosis of the respective rate 
distribution. 

Finally, Figures 5D-F demonstrate clearly how pattern for- 
mation occurs for intermediate /-values /™ d < /J nter < f& if the 
system is neither clearly in the mean- nor strongly fluctuation- 
driven regime, and even changes from one to the other regime 
with increasing coupling strength, cf. Figure 2C. As we demon- 
strate in the supplementary material section SI, most input and 
network settings have pattern formation onsets that usually agree 
much better with the /™ d prediction than with J™. This is a very 
consistent finding, even in cases where neurons are well in the 
fluctuation-driven regime in terms of subthreshold mean and 
pronounced input variance. 

As we discuss in the supplementary material section SI the rea- 
son for this finding lies in an asymmetry between the excitatory 
and inhibitory compound input current statistics. The excitatory 
subpopulation tends to be synchronized, even in the presence of 
strong balanced external noise, while inhibition actively decor- 
relates itself. The explanation for this effect lies in the local 
connectivity of ring networks, leading to a population-specific 
pattern of correlations already below the critical coupling: excess 
synchrony of the excitatory population will reinforce further spik- 
ing, while spikes from the inhibitory population decreases the 
instantaneous firing probability [see also Helias et al. (2013) for a 
related discussion of this effect in the framework of balanced ran- 
dom networks]. Near synchrony is thereby effectively enhanced 
in the excitatory population and suppressed in the inhibitory 
one. Volleys of excitatory input spikes act like compound pulses 
with large amplitude (on the order of up to 9) in an otherwise 
balanced asynchronous background activity. At these amplitudes 
the effective gain Equation (7) derived from linear response the- 
ory becomes linear in / with a slope proportional to 1/0, i.e., 
the slope of the linear model Equation (11), see supplementary 
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FIGURE 5 | The histograms of rates, variance and kurtosis for various 
weights J in a network of N = 2500, k = 250, and g = 6 averaged 
over 10 trials. (A-C) Shows the mean-driven case (/ x = 2500 pA 
delivered as Poisson noise with external synaptic strength J x = 0.1 mV), 
(D-F) the fluctuation-driven case (Poisson noise adapted s. t. \i[RI] = 5 mV 
and a[ffl] = 60mV) and (G-l) an intermediate case, where the external 
input was Poisson noise (current-amplitude 875pA, J, = 0.1 mV). For 
subcritical J < J c the distributions are approximately Gaussian with small 
variance and kurtosis close to zero. For supracritical weight the 
distributions become broader and platykurtic, indicated by the negative 
kurtosis. The vertical dashed lines in (B,C) and (E,F) mark the respective 



critical weights, i.e., J c = 0.506 mV [red, using Equation (13)] for the 
mean-driven, and J c = 0.905 mV for the strongly fluctuation-driven case 
[black, using Equation (9)1. The dashed vertical lines in (H,l) mark the 
respective predictions for the critical weight by Equation (13) (red, 
independent of r|), and by Equation (9) with only the contribution of 
3v/3|i (gray) and both linear and quadratic contributions Bv/S\l, dv/Sa 
(black) for the estimation of the effective coupling strengths, cf. Equation 
(7b). The increase in the kurtosis and decrease in the variance for very 
large J in (B,C) and (H,l) is due to the rectification of rates that leads to 
increasing mass at zero. Note that sampling is denser around the 
expected J c -value in the first two rows. 



material section SI. This explains the fact that Equation (13) has 
better predictive power also in the intermediate or — comparably 
weakly — fluctuation-driven scenarios. Moreover, the fluctuation- 
driven prediction only becomes valid if the additional external 
noise sufficiently dilutes residual synchrony, such as it is the case 
for [i[RI] = 5 mV and o[RT] = 60 mV. 

3.4. COARSE-GRAINED RING NETWORK 

In section 3.2 we saw that the full system can effectively be reduced 
to a five-dimensional system. However, the computation of eigen- 
systems (cf. Appendix A3) is still quite involved. In this section we 



study how well a further simplified system that is formally identi- 
cal to the well-known Ermentrout-Cowan networks [Ermentrout 
and Cowan (1979a,b, 1980) and supplementary material section 
S2] predicts the dynamics of the full network. 

We coarse-grain the system and combine groups 
of I = N/Ni neurons, such that they contain four 
excitatory and one inhibitory neuron (I needs to 
be odd in the following), as indicated in Figure 6. 
The neurons i c € {0, . . . , I — 1} within each cell 
c G {0, . . . , C — 1}, C = N/d, are connected to every other 
neuron within their cell c but not to themselves. Moreover, all 
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FIGURE 6 | We coarse-grain the system by introduction of cells 
containing four excitatory neurons and one inhibitory neuron 
[indicated in (A) by the dashed line]. The coarse-grained system is 
invariant to three symmetry operations. (A) T{ is a translation of the cell by 
I = 5, (B) R is the mirror symmetry of two neighboring excitatory neurons 
within each cell, and (C) S is the mirror symmetry of the two pairs of 
neighboring excitatory neurons within each cell. 



neurons within one local cell c are connected to all neurons within 
the K < C/2 neighboring cells to the left and to the right, i.e., 
{(c — K) mod C, . . . , (c — 1) mod C, (c + 1) mod C, . . . , (c + 
K)modC}. The computation of the eigensystem is outlined in 
Appendix A4. Due to the two mirror symmetries corresponding 
to exchange of neurons within the unit cell (with corresponding 
eigenvalues r, s € {— 1, 1}) and the translational symmetry with 
the eigenvalue determined by the wavenumber a, the system can 
be reduced to effectively two dimensions. 

3.4. 1. Coarse-grained network dynamics 

In the coarse-grained system we only need knowledge about two 
elements of the ring, one excitatory and one inhibitory neuron 
within the same cell (arbitrarily chosen as vq, V(^_i)/2 in c = 0, 
without loss of generality). For a fixed set of eigenvalues r, s, a, 
the activities of all remaining neurons are then uniquely deter- 
mined. The homogeneous mode v 0 (f) couples to the symmetric 
subspace r, s = 1 only and we can write v 0 (t) thus as a linear com- 
bination of the two eigenvectors v* 0 and v\ x , where arj = 0 (cf. 
Appendix A4). 

To analyze the stability of this mode with respect to spatial 
perturbations we write the population activity as a linear combi- 
nation of eigenvectors with non-zero wavenumbers a;, , a ., vj,., 
such that 



v(t)= E Ki^+O 



(18) 



i = 0 



Then, for all neurons n e {0, . . . , N — 1} the input, (Wv(t)) [ri], 
in the reduced system becomes 

(J[n] + K)modC I- 1 

(Wv(t)) [n] = J2 E ( W n,(cemodN+j) v(f)) 

c=(J[n]-K)modC j = 0 



\c=-K , = 0 I 



[J\n\l modN] 

-sj{ E T iJ2 fl <«<) 

\c = -K i = 0 I 

[J[n\lmodN +{l - l)/2] 

/C- 1 K 



= 41 ( E E ^ fl < ^ V U,a,) m0dN l 
\i = 0c = -K I 

-^(e E eiaiC 

\i = 0c = -K / 

[J[n]lmodN+ (I - l)/2] 
= 47 ( E Q[a ' ] fl ' (f) v U,a^ UlnV modN] 

(E "M^W^J l^modN 
+ (£-l)/2] 

with J^fn] := \n/l\. In the first identity we made use of the 
operators R, S and 7^ to express the input-connectivity of the net- 
work. The second identity follows from the linearity of the sum 
over eigenmodes and the respective eigenvalues of the operators 
R, S and Ti. We excluded synaptic self-coupling, so we need to 
subtract that input of the cell to itself, i.e. 

7 ^E a >'( f ) V U,a;^M ' ^ ~ K,J[n]lmodN+(e- l)/2) 

• On,J r [nKmodN+(f- l)/2- 



Leaving out one inhibitory input if the receiving neuron is 
inhibitory, and one excitatory input if the receiver is excita- 
tory leads to an asymmetry in the total input of excitatory and 
inhibitory neurons. This implies that the spatially homogeneous 
mode is not an eigenmode of the system anymore. To compensate 
for this the weights are rescaled such that the input per neu- 
ron equals that in the full (non-coarse-grained) system, i.e., to 
u. = K/(P-g(l-P)). Hence, 



J((3 + 8K)-g{2K+l)) 



=: [L £ W ijli€£ (20) 



m+(4-g)2K) 



[J[n]£ mod N + j] 



(19) 



Note, that if we do not exclude synaptic self-coupling in the 
coarse-grained system the eigensystem looks dramatically differ- 
ent from that of the full system. In particular the eigenvalues 
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are nearly fully real and the critical eigenvalue is much smaller 
implying a higher critical weight. 

In the homogeneous mode the input to each neuron n e 
{0, . . . , N — 1} is completely identical, all neurons fire at the same 
rate v 0 given by the self-consistent solution 



tmV„[0] 

t m v 0 [2] 



-6 0 \ U\L £ n[a\- 
0-9/ \ 4\ix£l[a] 



-glii(Q[a] - 1) 



RL 



9/2 
■9/2 



where £2 (a) := Sln(a( f +\ /2)) . 

v ' sin(a/2) 

rescaled coupling matrix 



The eigenvalues of the reduced 



W J lA^ £ Q[a] - 1 
4[xi^[a] 



-gV.x(Q[a] - 1) 



(21) 



then determine for which parameters g, J the spatially homoge- 
neous state becomes unstable and which wavenumber A,- (a,) will 
grow fastest and define the spatial pattern. The critical eigenvalue 
for the coarse-grained system is thus given by the eigenvalue of 

^g- with largest real part that will cross unity from below first, i.e. 

K(g,J)= max f— — (|x I (4f2[oi,]-l)-^ £ (f2[ a ,]-l) 

is {o c-i) I 26|ii|j,£ V 

+ y4gn, x n, £ (l - 5Q[a,]) + (indict,-] - 1) - gu. £ (£2[a ; ] - 1)) 2 )}(22) 

In the fluctuation-driven regime / and g in Equation (21) need 
to be substituted by / := W(J) and g, such that the critical 
eigenvalue is given by X^ d (g, /). 

On a mesoscopic level the coarse-grained and full sys- 
tem give rise to very similar activity patterns — differences 
become apparent in the fine-scale features, e.g., the full 
analysis of a network of 2500 neurons with 10% con- 
nectivity and g = 6 predicts destabilization of wavenum- 
ber 14 in both mean- and fluctuation-driven regime, while 
the coarse-grained model predicts 13, see also Figure 7. 
We note that the coarse-grained model presented in this section 
is formally identical to the linearization derived in Ermentrout 
and Cowan (1980) in the context of a neural field ring model 
if one assumes a boxcar-footprint. We summarize the respective 
derivation in the supplementary material section S2. 

3.5. SENSITIVITY TO RANDOMNESS 

To conclude, we study how robust the findings of the previ- 
ous sections are to effects of randomness in either structure or 
weight distribution. Though we again only show the eigenvalue 
spectra for the synaptic coupling matrix W/9, that is the rel- 
evant linear operator in the mean-driven case, all results map 
straight-forwardly to the fluctuation-driven regime with appro- 
priate rescaling of / and g. 

3.5. 1. Small-world networks: structural noise 

First, we will study the effect of the introduction of structural 
randomness by rewiring individual connections with proba- 
bility p r [for details see Kriener et al. (2009)] such that the 
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Wavenumber A =No/2ir£ 

FIGURE 7 | (A) Sketch of the layout of the coarse-grained system. Each 
neuron is connected to the L — 1 other cells within its cell and to all neurons 
in its 2K nearest cell neighbors. A shift in neuron label by five yields the 
same ring as before and is formally expressed by the shift operator Ji with 
1 = 5. (B) Shows the coupling matrix W of the ring: black squares depict 
inhibitory, white squares excitatory, and gray squares zero coupling from / 
to / (W = 60, k = 30). (C) Eigenvalue spectrum of a scaled ring matrix W/8 
with N = 2500, k = 250, and g = 6. (D) The two eigenvectors belonging to 
the twice degenerated eigenvalue from (C) with largest real part, XJT d . (E) 
Five bands of eigenvalues as a function of the wavenumber. The black line 
depicts the band containing X md . The inset shows the critical bands of the 
system in section 3.2 (black circles, cf. also Figure 3) vs. the critical band in 
the coarse-grained system (gray triangles). The critical wavenumber differs 
by one. (F) The corresponding rate per neuron in a simulation of W = 2500 
integrate-and-fire neurons with J = 1 mV and l x = 750 pA. 



input composition from the network to every neuron stays con- 
stant at u,. The random rewiring procedure means that each 
realization of a network will be different and moreover lacks 
invariance to Tg such that the eigensystems of individual cou- 
pling matrices need to be computed numerically. For low p r 
the rate model for the translation-invariant network still gives 
very good results and the patterns that form for large / are 
indeed strongly correlated with the critical eigenvector v c , see 
Figure 8. For higher p r predictions become worse and more than 
one vector can contribute to the pattern (e.g., Figures 8C2,C3). 
For low rewiring probability p r < 0.02 we can still analyt- 
ically estimate the critical eigenvalue by employing a mean 
field model (Grabow et al., 2012) for small-world networks 
in which we rewire "on average," cf. Appendix A5 and 
Figures 8A2-C2,A4-C4. Hence, also for small-world networks 
of excitatory and inhibitory spiking neurons, rates and linear 
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FIGURE 8 I Effect of rewiring of edges on the eigenvalue distributions 
(row 1), the mean-field density of the real parts of eigenvalues (row 2), 
the crictical eigenvectors of the specific realization (row 3) and the 
mean-field system (row 4) and the rate per neuron (row 5) in 
small-world networks for three different rewiring probabilities 
p r = 0.01, 0.03, 0.05 in (A,B,C), respectively (cf. Watts and Strogatz, 
1998; Kriener et al., 2009). While in (A3.B3) it is the dominant eigenmode v c 
that correlates strongest with the actual rate distribution [(A5.B5), Pearson 



correlation coefficients of c = 0.967 and c = 0.945, respectively], in (C3) the 
linear combination of the two leading eigenvectors correlates best with the 
actual rates l(C5), v c has c = 0.68 while the combination of the two leading 
modes gives c = 0.91 ]. The inset in (A2) shows the estimate for the critical 
weight J c as a function of p r in the mean field model (MF, black) and as an 
average of 50 network realizations (SW, gray). Parameters: N = 2500, 
J = 1 mV and l x = 750 pA. During rewiring the number of excitatory and 
inhibitory inputs was kept constant. 
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stability can be estimated from the translation-invariant network 
rate model if p r is small. The results of the mean-field model are 
slightly worse than for unweighted networks (cf. Grabow et al., 
2012) because of the two qualitatively different types of edges 
in the networks considered here (excitatory and inhibitory). For 
the small number of edges k we assume in Figure 8 it can occur 
that for very small p r no inhibitory edges are rewired at all, while 
the mean-field model still reassigns inhibitory connection den- 
sity. We find the mean field to work better for denser or larger 
networks, and if the excitatory and inhibitory input per neuron is 
kept the same during rewiring (not shown). 

If the constant input per neuron constraint is dropped, the 
homogeneous mode (1, 1) T is not an eigenmode anymore and 
also the prediction of stability fails in large parts. On the one hand 
there are activity inhomogeneities forming out simply due to the 
fact that some areas on the ring will by chance get more or less net 
inhibitory input, and hence some neuronal subpopulation will be 
driven to zero rates, while others — lacking net inhibitory input 
from these rectified neurons — will have very high rates. Similar 
effects arise when the distribution of inhibitory neurons is not 
periodic anymore, but, e.g., distributed randomly along the ring. 

3.5.2. Ring networks not conform with Dale's principle 

Finally, if Dale's principle — i.e., the biological fact that each neu- 
ron can only either depolarize or hyperpolarize all its postsynaptic 
targets, but never both at the same time — is violated, both the 
eigensystem (Figure 9A), as well as the spiking dynamics looks 
akin to that of a random network, even if the underlying topol- 
ogy is a ring, cf. Figures 9B,C. Hence, the identity of neurons in 
terms of being excitatory or inhibitory is a necessary condition 
for the formation of structured patterns in the case of identical 
footprints for inhibition and excitation considered here. 

Even though in all cases pattern formation is either strongly or 
completely impaired, the prediction Equation (14) of individual 
rates in the low coupling regime is excellent, also in inhomoge- 
neous networks, where different neurons receive very different 
inputs, given the individual mean membrane potentials (V,(f)) f 
are used in Equation (15). 



4. DISCUSSION 

In this paper we analyzed the dynamics of pattern formation 
in networks of excitatory and inhibitory spiking neurons that 
are embedded on a ring architecture. In particular, we studied 
the dependence of the underlying bifurcation on the input cur- 
rent statistics and structural noise, and derived a coarse-grained 
system that is formally analogous to the classical Ermentrout- 
Cowan networks. To conclude, we will summarize and discuss our 
findings. 

4.1. IMPACT OF INPUT CURRENT STATISTICS 

The very nature of spiking neuron networks leads to consider- 
able input current fluctuations, and hence in general both mean 
|x[_Rf] and standard deviation o[RI] of the input will shape the 
output rate of the neuron. For leaky integrate-and-fire neurons 
the output rate as a function of [i[RI] and o[RI], assuming sta- 
tionary Poisson-like uncorrelated input spike trains and weak 
synaptic coupling /, is given by the Siegert equation (5) (Amit and 
Tsodyks, 1991; Amit and Brunei, 1997). Thus, in order to analyze 
the linear stability of this self-consistent solution, in general the 
derivatives with respect to both [i[RI] and o[RI] need to be taken 
into account, cf. Equation (7). 

We analyzed two different ways to drive the networks: in the 
first all neurons in the network were driven with the same exci- 
tatory Poisson-noise vx, independent on the recurrent coupling 
strength /. In this case the total input RI the neurons receive, i.e., 
the external input RI X together with the recurrent input _R7 S from 
the network, changes with increasing /, such that the mean \l[RI] 
decreases, while the standard deviation a[RI] increases. 

In the second input scenario, we corrected for the change 
in recurrent input contributions with changing / by adminis- 
tering compensating inhibitory and excitatory external Poisson 
input vix, ve x . Thus, \l[RI] and a[RI], and hence also v 0 are 
approximately constant with increasing /. 

We find that in the first input scenario the critical J c as 
observed in simulations deviates strongly from the one predicted 
by Equations (7), (9), which are derived from the self-consistent 
solution Equation (5) by linear response theory. In particular, 
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FIGURE 9 | (A) Eigenvalue distribution in the complex plane of a ring network 
with random assignment of weights, irrespective of the identity of a neuron, 
hence each existing synapse has weight J with probability p and — gJ with 
probability (1 — P). These ring coupling matrices have eigenvalue spectra very 
akin to those of corresponding random networks and most of the 
eigenvalues adhere to the circle law prediction (Girko, 1984) for random 
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networks r = ^ Wp(p + gr 2 (1 — P)) — \l 2 /N (dark line). Only a few singular 
eigenvalues on the left of the center indicate the underlying ring topology. 
The eigenvectors of such "hybrid" ring matrices have no apparent structure 
and no pattern formation occurs when the system becomes supracritical, as 
clearly visible from the spiking activity (B) and the rates per neuron (C). Here, 
J c w 0.44 mV, simulation parameters J = 0.6, g = 6 and N = 2500. 
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for the mean-driven regime, where [i[RI] > V^, J™ becomes 
basically independent of the rate v 0 as well as of a[RI], irre- 
spective of its yet considerable magnitude. Instead, a very simple 
linearization derived from the noiseless /-/-curve Equation (10) 
applies. The /™ d derived from this linearization is always con- 
siderably smaller than the prediction by Equations (7), (9). For 
intermediate cases where neurons are neither mean- nor strongly 
fluctuation-driven, the critical coupling strength lies in between 
the two predictions, and often closer to /™ d . 

Possible reasons for the failure of the prediction of Equations 
(7), (9) in the constant external drive scenario (parameterized 
by T), cf. section 2) are the break-down of the weak-coupling 
assumption, as well as the uncorrelated Poisson assumption for 
the input spike trains in the recurrent contribution. In line with 
earlier findings (Tetzlaff et al., 2012), we could not identify devi- 
ations with respect to coupling strength, however, we observed 
clear deviations from the Poisson assumption. 

If spiking is Poissonian we expect the interspike intervals to 
be distributed exponentially, and the coefficient of variation (CV) 
of the interspike intervals (ISI) to be unity. A smaller CV indi- 
cates more regular, a larger one more irregular spiking. As can be 
seen from Figures 1A-C and 1G-I, activity is spatio-temporally 
structured, with locally clustered spiking that can induce signifi- 
cant pairwise correlations, and low coefficients of variation of the 
interspike intervals for small / (cf. also supplementary material 
section SI, Figures S1B,C). In effect, these deviations from uncor- 
related Poisson spiking lead to deviations from the Gaussian white 
noise approximation underlying the derivation of Equations (5) 
and (9), especially for increasing /. 

We observe that in particular the excitatory subpopulation 
tends to synchronize because of the highly recurrent local struc- 
ture and positive feedback, while the inhibitory population 
actively desynchronizes itself due to negative feedback (see sup- 
plementary material section SI, Figures S3, S4, and Helias et al. 
(2013) for a related discussion). Thus, high amplitude volleys of 
excitatory input in a balanced background comprised of asyn- 
chronous inhibition and external Poisson-drive move the net- 
work in a regime that is better captured by the mean-driven 
low-noise approximation Equation (11). 

Activity in the strongly fluctuation-driven regime, on the other 
hand, is much more asynchronous and irregular by construc- 
tion, cf. Figures 1D-F with CV[ISI] of unity and larger (see 
supplementary material, section SI, Figures S1B,C). When in this 
scenario a[RI] is very large, while \n[RI] is sufficiently subthresh- 
old, the prediction for the critical coupling strength ]~ from 
Equations (7), (9) for linearity to break down indeed appears to 
be correct, or even seems to underestimate the actual onset of 
pattern formation. 

This shallow transition around can be explained by slow 
noise-induced transitions between the two degenerate eigenvec- 
tors which grow quickest once pattern formation takes place. 
However, the fact that the CV is larger than unity will also lead to 
deviations from the prediction. In a series of papers Moreno et al. 
(2002); Renart et al. (2007) and Moreno-Bote et al. (2008) stud- 
ied the impact of deviations of the Poisson-assumption of input 
spike trains on the firing rates of neurons. They studied the effect 
of positive and negative spike train auto- and cross-correlations 



parametrized by the Fano-factor F = a 2 [counts] /u, [counts] of 
spike counts on the output firing rate of current-based LIF neu- 
rons in both the fluctuation- and the mean-driven regime. For 
renewal processes the CV of interspike intervals is directly related 
to the Fano-factor of counts by _F = CV 2 (Cox, 1962). In Moreno 
et al. (2002); Renart et al. (2007); Moreno-Bote et al. (2008) it 
was shown that in the fluctuation-driven regime output rates are 
quite sensitive to input noise correlations (parameterized by the 
Fano factor), while in the mean-driven regime sensitivity is less 
pronounced. 

In simulations of ring networks in the strongly fluctuation- 
driven case we observe that rates decrease with increasing cor- 
relations, while in otherwise equivalent random networks, rates 
increase (not shown). So the different recurrent input structure 
plays a crucial role for the effective input current and thus net- 
work dynamics properties. These differences are often, as here for 
ring and grid networks, a direct consequence of complex network 
topology, and it is thus important to understand how connectivity 
structure translates to single neuron activity, as well as collective 
network states. A more thorough analysis of such effects in the 
context of firing rate stability in random and spatially structured 
networks will be subject of subsequent research. 

4.2. IMPACT OF FINITE TIME SCALES 

We assumed throughout the manuscript that the system has a 
negligible refractory time constant, small delays and instanta- 
neous post-synaptic currents (8-synapses). If these assumptions 
are dropped we observe differences especially for the mean-driven 
case: for larger delays d in the order of several milli-seconds, 
pronounced oscillations can occur that tend to synchronize the 
system locally on the order of the footprint (Kriener et al, 2009), 
or for very strong external drive also on a population-wide level 
(Brunei, 2000). A systematic study of the complex effects of delays 
on pattern formation in neural-field-type ring networks of spik- 
ing neurons was presented in Roxin et al. (2005). In the network 
we analyzed here, noisy delay oscillations mainly perturb develop- 
ing patterns and can thus shift the onset of clear pattern formation 
to larger J c . 

Similar effects arise in the case of finite refractory time T. r ef- 
A larger i: re f, especially in the presence of oscillations, i.e., highly 
coincidental input spiking, decreases the effective network input, 
in particular if t re f > d, and thus also shifts the input regime. This 
can be compensated for by explicitly taking into account x re f in 
the derivation of the noiseless linearization. 

Finally, if synaptic time constants are finite, while total input 
currents stay the same, delay oscillations can be smoothened out, 
and the effect of absolute refractoriness is decreased. Increasing 
synaptic time constants can thus counteract the effect of finite T re f 
and d. 

4.3. SPIKING NETWORKS vs. NEURAL FIELD MODELS 

Pattern formation, such as activity bump formation or peri- 
odic patterns, is a well-studied phenomenon in ring and toroidal 
networks (see e.g., Ermentrout and Cowan, 1979a,b; Ben-Yishai 
et al., 1995; Usher et al., 1995; Bressloff and Coombes, 1998, 2000; 
Shriki et al, 2003; Marti and Rinzel, 2013). Earlier studies of 
periodic pattern formation in spiking neuron networks already 
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showed that ring- and grid-networks of mean-driven neurons can 
undergo a Turing-bifurcation (Usher et al., 1995; Bressloff and 
Coombes, 1998, 2000). In particular, Usher et al. (1995) observed 
a dependence on the input level. For weaker, yet suprathreshold 
mean input current, activity patches diffuse chaotically across the 
spatial extend of the system, while for stronger inputs a stable 
spatial pattern develops that relates to the critical eigenmode. We 
observe similar dynamics in the networks considered here as well, 
if the system is in the intermediate regime. This can be understood 
as follows: if the noise component is relatively strong, it shifts 
the actual critical coupling strength /™ d to higher values /' nter 
inbetween /™ d and J™. If the system is thus effectively slightly 
sub-critical in terms of /™ ter , input-noise can transiently excite 
eigenmodes close to effective criticality, due to Hebbian (Dayan 
and Abbott, 2001) or non-normal amplification (Murphy and 
Miller, 2009). 

To compare our results to the classical work on neural field 
models on ring topologies (Ermentrout and Cowan, 1979a,b, 
1980) we moreover studied a coarse-grained model that reduces 
the AT-dimensional system to an effectively two-dimensional one 
that is structurally equivalent to the linearized two-population 
model presented, e.g., in (Ermentrout and Cowan, 1980) (see 
supplementary material section S2) which, however, allows for a 
direct quantitative mapping to the spiking network dynamics. 

Because of the coarse-graining some of the connectivity details 
get lost and thus the critical eigenmode is not identical with that 
of the original ring network. Moreover, we observe the impor- 
tance of excluding self-coupling in the coarse-grained model: if 
that is not excluded the resulting eigensystem looks dramatically 
different with nearly completely real-valued eigenvalues. 

4.4. IMPACT OF STRUCTURAL NOISE 

For the translation-invariant distribution of inhibitory neurons 
across a ring topology the eigensystem, and hence the respective 
linear stability, can be computed analytically. If inhibitory neu- 
rons are however randomly distributed, some parts of the ring will 
have a higher, others a lower density of inhibitory neurons, such 
that these networks form an ensemble of possible realizations, 
most of which are not invariant to translations. The resulting 
inhomogeneous rate per neuron in the subcritical regime can still 
be predicted well from the linear rate model. However, even when 
knowing the dominant eigenvector, it will usually not be peri- 
odic, and due to the inhomogeneous input per neuron different 
neurons will be at quite different working points. There is inter- 
ference between the rate distribution and the dominant mode. 
Thus, the activity pattern that evolves in the supracritical regime 
is not straightforward to predict. 

If the weights are distributed fully randomly across the whole 
ring topology irrespective of the identities of the neurons, i.e., in 
violation of Dale's principle, the resulting eigensystems look very 
similar to those of corresponding networks with random topolo- 
gies. Only a few surviving singular eigenvalues serve to distinguish 
the underlying ring from the random topology. This directly 
translates to the spiking dynamics: the activity is much more asyn- 
chronous than that of Dale-conform networks and most of all 
pattern formation can not take place due to the effective lack of 
lateral inhibition. 



4.5. IMPACT OF NETWORK SIZE 

If k/N is kept constant when varying the network size N the 
critical coupling strength in both the mean- and the fluctuation- 
driven regime decreases. For example, for N = 10000 and k = 
1000 the critical coupling strength becomes /J™ 1 » 0.2 mV for 
the mean-driven, and ]^ «a 0.32 mV for the fluctuation-driven 
regime. If the number of synapses per neuron k is fixed when N 
increases, both /™ d and J™ hardly change because the maximal 
eigenvalue is approximately unaffected by network dilution. 

We emphasize that the analysis presented here can be extended 
to two-dimensional grid networks in a straightforward manner, 
see supplementary material section S3. We conclude by remarking 
that the rate model Equation (13) is also a valuable tool to study 
the rate dynamics of spatially embedded spiking neuron networks 
with e.g., inhomogeneous neuron distribution and distance- 
dependent connectivity in the mean-driven regime, allowing for a 
treatment of more general networks than commonly studied ran- 
dom networks or neural field models. The presented model is in 
conclusion a further step in the efforts to find and understand 
appropriate descriptions of the high-dimensional spiking activity 
in structured networks. 
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A. APPENDIX 

A1. MODEL AND SIMULATION DESCRIPTION 

Table A1 | Model description, simulation paradigm and parameters. 







Populations 


Three: excitatory, inhibitory, external input 


Connectivity 


Ring connectivity with excitatory and inhibitory neurons arranged on the same ring, such that each fifth neuron is inhibitory; 
same connection footprint for both populations 


Neuron model 


Leaky integrate-and-fire (LIF), fixed voltage threshold 
precise spike timing (Hanuschkin et al., 2010) 


Synapse model 


8-current inputs (discontinuous voltage jumps) 


Input 


Independent Poisson spike trains 






Name 


Elements Size 


E 


LIF neuron W E = 4A/| 


I LIF neuron Ni 


E x 


excitatory Poisson input generator one realization per neuron 


lx 


inhibitory Poisson input generator (only for the fluctuation-driven case) one realization per neuron 




C 


CONNECTIVITY 


Name 


Source Target Pattern 


EE 


E E connect to all excitatory neurons within a given distance on the ring, weight J, delay d 


IE 


E I connect to all excitatory neurons within a given distance on the ring, weight J, delay d 


El 


I E connect to all excitatory neurons within a given distance on the ring, weight — gJ, delay d 


II 


I I connect to all excitatory neurons within a given distance on the ring, weight — gJ, delay d 


E x 


E x E U I independent Poisson spike trains, weight J x , rate vg x 


lx 


l x E U I independent Poisson spike trains, weight —gJx, rate V| x 



(Continued) 
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Table A1 | Continued 







Name 


LIF neuron 






Type 


Leaky integrate-and-fire, 8-current input 






Subthreshold dynamics 


Xm V(t) = -V(t) + Rl(t) 

V(t) = Vres 

\A/+A/ ext x — , 

R/(f) = E /=1 w ^ w a T, k - - d » + R, * ( f) ■ 


if 

else 

Wij 


t > t* + x ref 
6 {-gJ, 0, J] 


Spiking 


if V(r-) < e a V(t+) > e 

1. set f* = f 

2. emit spike with time-stamp t* 







Type 

Poisson generators 



INPUT 

Description 



Rate vg x , v| x , each generator projects independent realizations to all neurons 



A2. APPROXIMATION OF THE SIEGERT FORMULA BY AN 
AFFINE-LINEAR INPUT-OUTPUT RELATION 

For strong mean drive \i[RI] S> 6 we can approximate the input- 
output relation ( 10) by an affine linear relationship. In the follow- 
ing we assume the firing rate 0 < v <SC 1 /t re f, which is the case for 
all relevant regimes, so we can neglect the effect of the refractory 
time t re f. We use the Ansatz 

v(D = - (t m Log [1 - e/(itf)]) -1 =a+bRI (Al) 

<3> -t m Log [1 - 9/(R7)] = (a + bRI)' 1 

By taking the derivative with respect to RI on both sides of the last 
expression we obtain 



1 - Q/(RI) 



Q/(RI) 2 = -b(a + bRiy 



(A2) 



6t„ 



QRI-(RI) 2 a 2 /b + 2aRI + b(RI) 2 



Equating the terms in proportion to (RI) and in proportion to 



(RI) in the denominator fixes the parameters 

l 



b = ^- , leading to the approximation 



2t n 



v(RI) 



1 1 
+ RI 

2T m Otrvi 



Indeed, to show that for 9 » RI 
1 



tniLoi 



1 1 

H RI. 



and 



(A3) 



(A4) 



we introduce the auxiliary variable Q/RI =: (1 — x) and compute 
the limit 



lim 

x/ 1 



1 



+ 



1 



LogM 1 



: lim 

%/ 1 



1 — x + LogM 
Log[x](l -x) 



(AS) 



Both numerator and denominator of this expression tend to zero, 
thus we apply L'Hopital's rule and obtain for (A5) 



lim 

x/ 1 



I/X- 1 



1 jx — 1 — Log[x] 



(A6) 



Again, numerator and denominator tend to zero, thus we apply 
L'Hopital's rule again and obtain for (A6) 



lim 

x/l 



l/x 2 



(l+x)/x 2 _ 



■■ lim 

x/\ 



1 +X 



(A7) 



This yields (A4). 

As a heuristic observation, we note moreover, that in the large- 
r-limit , i.e., a 0 ^> [L 0 , the integrand of Equation (5) becomes 
basically unity, and the integral thus scales as v^ 1 ~ *Jiix m Q/a 0 . 
In particular we observe that the input-output-relation Equation 
(5) for Vres = 0 takes the following linear form in a 0 : 



v 0 (a 0 , |x 0 ) 



tmv^ 



0o + M-o 



(A8) 



The input-output-relation of the LIF is thus also linear in a 0 for 
Oo S> M-o. see supplementary material section S4. 

A3. EIGENSPECTRUM OF THE RING GRAPH WITH EQUIDISTANT 
INHIBITORY NEURONS 

We consider graphs, in which each neuron i is connected to 
neurons i—\c/2 to i + k/2, k even, but not to itself. To keep 
analysis simple, we assume that each l-th neuron, N mod I = 0 
is inhibitory. 

Define Ti as the translation operator that shifts all neuron indices 
by I such that for every vector v we have 



Tt v = Te(vo, vi,v 2 , vn-i) 

= (Vl, V(+l, V e+2 , Vjv-l, v 0 , 



, n-i) T (A9) 
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i.e., for the node index n Injectivity: 

Ker(P ; or|) = {0} 

T t [ri\ = (n + l) modN , n e {0, N - 1} . (A10) 

Subjectivity: 

Ti is a unitary operator and has eigenvalues cp; and eigenvectors 

w h such that dim [(P, o r\)(C e )] = dim [Eig(T € , (p ; )] = I . 

Tgwi = e ,a, wi =: cp/w; , (All) 

Now, we can calculate the spectrum of W. Due to [ W, Ti ] = 0, i.e., 
w ith W and Ti commute, we can diagonalize both matrices, such that 

N/ . they have a common system of eigenvectors. Hence, it is sufficient 

T e wi = wi =>■ ct; = 2nU/N, I e {0, N/€ - 1} . (A12) to reduce the problem to the ^-dimensional vector-space embed- 
ded in C N , spanned by eigenvectors ofEig(Ze, e m )\ c t, i.e., ivj, Z € 
Hence the N/l Mold degenerated eigenvalues of T e are ip; = rn, - 1}. To obtain the eigenspaces Eig(W,X ;i ) we only 

e «2H«/N ; e {0; N /i _ i}. The eigenvectors w, are discrete need t0 solve the diagonalization problem for the corresponding 
harmonics and fullfil x £)- ma trices rT 1 o W o P, o r\ = (Pi W)\ c i*t : 

W W)modN = <P' W 'imodN ' ( A13 ) O W O P, O T))W = Xw (A20) 

We define Once we got the I eigenvalues and eigenvectors {X; ;, iv/.,} of W 

t j in the respective sub-spaces, the operator P; o t) generates the 

C 3w=(%wi,w 2 ,w 3 ,...,wn) (A14) correS ponding N-dimensional eigenvector w,. 

and the map ti that embeds the C l as a subspace into the C N A4 EIGENSYSTEM OF THE REDUCED RING SYSTEM 

We now consider a coarse-srained network model where con- 

i-\ 

£ l c f £ N ~ ~ ~ (A15) nectivity i s defined between C = N/£ cells of I neurons each as 

^ " C ' W ^ W ' e ' ~~ ' V ' depicted in Figure 6A. Each neuron is connected to all other I— 1 

neurons within the cell it resides in, as well as to all neurons in 
where e„ i e {0, N - 1} are the N canonical basis-vectors the neighboring 2K = pC cells, yielding block-wise connectivity 
of the R N and w, € C is the i-th vector component of w. We matnces as indicated in Fi g ure 7B - We introduce a new global 
moreover define the operator family P,, I e {0, N/l - 1} neuron labelm S such that 

N/t-i n = c£ + i c , c=\-\, and i c = nmodl (A21) 

P,:C f cC N ^C N , Y, T{e- i2 *W N v. (A16) UJ 

3 = 0 Symmetries: 

........ , . The system in invariant to shifts about elementary cells, i.e., for 

We see that (Pi o ri) indeed projects into the eigenspaces and is an ,, . , , . . . ... 

, ° the cell index n we have again invariance to the operation 

isomorphism 



P,or| :C 1 ^ Eig(7>,cp ; ) CC W (A17) 



Ti[n] = (n + £) mod N . ( A22) 



, ... . , , Moreover, there are two mirror symmetries R and S within each 

as can be readily checked: „ . . _. , , ' 

... . . , . , .... cell (cf. Figures 6B,C), such that 

We see that for each I e {0, N/l - 1} 6 



N/l-\ N/e-i 
TdPiori) = T e J2 T{e- i2 *M N = £ if 1 e~ i2 *W N R W[n] = Rw [(ct+k)modN] 

j=0 j=0 

N/e-i 



W[(rf+(i-; c ))modjV] if! c <l 

(A23) 

w [n] if i c > 1 



(*) e i2^/N T k ie -' 2 ^ N = ^,0^, (A18) 

*: = 0 SW[„] = SW[( rf+ ; c ) mo dN] = W[(ct+e-l-i c )modN] (A24) 

where in the equation marked by (*) we performed an index shift R and S are operators, such that for the eigenvalues rj and sj, 
and used (e- ,a Ti) N l l = (e-' a Ti)°. and for the eigenvectors Vj and Wj we have R 2 Vj = r 2 Vj = vj and 
Linearit y ; S 2 w, = sjwj = Wj. Hence, rj,Sj € {—1, 1}. We have three corn- 
er-, ,t, i ,t. x muting operators Tp, R, S, i.e., 
V v,w€C«, fle c (Pi or\) (a w) = a (Pi or\)w and (P;ot))(w + v) 

= (P, o n)w + (P/ o n) v (A19) [T/.R] = [T/.S] = [S,P] =0. 
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So, we can find a common basis of eigenvectors to span the com- 
plete system. For each wavenumber A,- = a,N/27t£ one eigenvec- 
tor is obviously given by 



possible eigenvalues of _R, S each yields 



N/l-\ 



J2 (e~ ia 'TtYv { i-i)i2e ( i-i)/2 (A25) 



where e(i-i)n is the (I — l)/2-th canonical basisvector of the R N 
and (I — l)/2 is the index of the first inhibitory neuron on the 
ring enumerated as in Figure 6, since this dimension does not 
couple to either R or S. To get the residual N — N\ eigenvectors, 
we choose 



c- 1 

^s,a, ■= J2 ( £, ~' a ' T » C ( 1 + sS )( J + rK ^ e o (A26) 

c=0 

with first canonical basis vector en- 

Indeed, we see again that for each c e {0, N/l — 1} 



Tde-^Ttf = e- ,ac T c + l = e ' a e-' a(c+1) r'+ l 
= e' a (e- ia T e ) c+1 
and (e- ia T e f = 1 , 



Tp- 1 (A27) 



and 



X (1 + xX) = X + xX 2 = X + xl = x(xX + 1) 

= x(l+xX) (A28) 



with x € {r, s} and X e {-R, S}, so V r ,s,a are eigenvectors with 
eigenvalues r, s, a. Evaluation of the eigenvectors for the two 



c- 1 



"l,l,a 



J2 {e-^TtYiivo, v 0 , 0, v 0 , vo, 0 0) T ) 

c = 0 
C- 1 

vE i.i.o, = J2 ( e_<a ^) C (( v 0, -vo, 0, -v 0 , vo, 0, . . . , 0) T ) 

c = 0 
C- 1 

£ (e-' a ^) c ((v 0 , v 0 , 0, -v 0 , -vo, 0, . . . , 0) T ) 

c = 0 
C- 1 

^ (e-^T^^Cvb, -v 0 , 0, v 0 , -vo, 0, . . . , 0) T )(A29) 



"l.-l.o, 



-1,-1, a 



c = 0 



A5. MEAN-FIELD APPROXIMATION OF SMALL-WORLD NETWORKS OF 
EXCITATORY AND INHIBITORY NEURONS 

Analogous to the mean-field approach developed in (Grabow 
et al., 2012) for unweighted networks we will now compute the 
eigensystem of the mean-field small world network derived from 
the ring network with equidistant inhibition by random rewiring 
of edges (Kriener et al., 2009). In this case instead of assigning 
explicit weights W» G {0, /, — gj} to individual entries of the cou- 
pling matrix, all possible connections carry their expected weight 
as a function of the rewiring probability p r such that all weights 
within the original ring neighborhood are scaled by pi (p r ), and all 
weights that were 0 in the original ring are set to Wy = p2 (p r ) J if 
j is excitatory, and W y = —pi(p r )g! if j is inhibitory, where (cf. 
Kriener et al, 2009; Grabow et al, 2012) 



Pi(pr) = (1 ~Pr) + 



P 2 K 



(N-I-Kd-Pr)) 



(A30) 



Plipr) 



(N-l-K(l-p r )) 



The resulting average coupling matrix is again invariant with 
respect to Ti and the eigensystem can be determined exactly as 
described in Appendix A3. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



January 2014 | Volume 7 | Article 187 | 21 



